In [7]:
%pip install -q --no-warn-conflicts malariagen_data
import warnings
warnings.simplefilter(action='ignore', category=Warning)
import malariagen_data
import numpy as np
import pandas as pd
import os

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 165.1/165.1 kB 2.8 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 3.2/3.2 MB 38.0 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 35.8/35.8 MB 32.0 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.5/7.5 MB 82.9 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.0/4.0 MB 84.4 MB/s eta 0:00:00
  Preparing metadata (setup.py) ... done
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 302.5/302.5 kB 29.5 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 141.1/141.1 kB 14.5 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 24.4/24.4 MB 17.3 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.1/8.1 MB 61.7 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 210.2/210.2 kB 21.8 MB/s eta 0:00:00
  Preparing metadata (setup.py) ... done
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.6/8.6 MB 62.3 MB/s eta 0:00:00
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.6/1.6 MB 71.5 MB/s eta 0:00:00
  Building wheel for dash-cytoscape (setup.py) ... done
  Building wheel for asciitree (setup.py) ... done

mount google drive

In [ ]:
!pip install nbconvert
In [8]:
try:
    from google.colab import drive
    drive.mount("drive")
except ImportError:
    pass

Drive already mounted at drive; to attempt to forcibly remount, call drive.mount("drive", force_remount=True).

define directory within google drive as data cache

In [9]:
results_dir = "drive/MyDrive/malariagen_data_cache"
In [10]:
ag3 = malariagen_data.Ag3(results_cache=results_dir)
ag3
Out[10]:
MalariaGEN Ag3 API client
Please note that data are subject to terms of use, for more information see the MalariaGEN website or contact support@malariagen.net. See also the Ag3 API docs.
Storage URL gs://vo_agam_release/
Data releases available 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 3.10
Results cache /content/drive/MyDrive/malariagen_data_cache
Cohorts analysis 20240418
AIM analysis 20220528
Site filters analysis dt_20200416
Software version malariagen_data 11.0.0
Client location Utah, United States (Google Cloud us-west3)

define sample sets

In [11]:
sample_sets = ["AG1000G-BF-A", "AG1000G-BF-B", "AG1000G-BF-C"]

GAMBIAE¶

2L¶

In [12]:
sample_query="cohort_admin2_year == ['BF-09_Houet_gamb_2014', 'BF-09_Houet_gamb_2012']"
contig="2L"
In [13]:
ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)
In [14]:
window=2500
In [15]:
ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)










    
    


















    
    







  











    
    
































In [16]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_gambiae_Chr2L.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_gambiae_Chr2L.tsv')

# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_gambiae_Chr2L.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_gambiae_Chr2L.tsv




























GAMBIAE¶
2R¶
figure out window size



















In [17]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_gamb_2014', 'BF-09_Houet_gamb_2012']"
contig="2R"



     



















In [18]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [19]:



     


window=1000



     



















In [20]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [21]:



     


# Parameters

output_file_path = os.path.join(results_dir, 'H12counts_gambiae_Chr2R.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_gambiae_Chr2R.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_gambiae_Chr2R.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_gambiae_Chr2R.tsv




























GAMBIAE¶
3L¶


















In [22]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_gamb_2014', 'BF-09_Houet_gamb_2012']"
contig="3L"



     



















In [23]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [24]:



     


window=1000



     



















In [25]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [26]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_gambiae_Chr3L.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_gambiae_Chr3L.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_gambiae_Chr3L.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_gambiae_Chr3L.tsv




























GAMBIAE¶
3R¶


















In [27]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_gamb_2014', 'BF-09_Houet_gamb_2012']"
contig="3R"



     



















In [28]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [29]:



     


window=1000



     



















In [30]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [31]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_gambiae_Chr3R.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_gambiae_Chr3R.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_gambiae_Chr3R.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_gambiae_Chr3R.tsv




























GAMBIAE¶
X¶


















In [32]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_gamb_2014', 'BF-09_Houet_gamb_2012']"
contig="X"



     



















In [33]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [34]:



     


window=2500



     



















In [35]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [36]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_gambiae_ChrX.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_gambiae_ChrX.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_gambiae_ChrX.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_gambiae_ChrX.tsv




























Coluzzii¶
2L¶


















In [37]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_colu_2014', 'BF-09_Houet_colu_2012']"
contig="2L"



     



















In [38]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [39]:



     


window=1000



     



















In [40]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [41]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_coluzzii_Chr2L.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_coluzzii_Chr2L.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_coluzzii_Chr2L.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_coluzzii_Chr2L.tsv




























COLUZZII¶
2R¶


















In [42]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_colu_2014', 'BF-09_Houet_colu_2012']"
contig="2R"



     



















In [43]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [44]:



     


window=1000



     



















In [45]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [46]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_coluzzii_Chr2R.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_coluzzii_Chr2R.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_coluzzii_Chr2R.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_coluzzii_Chr2R.tsv




























COLUZZII¶
3L¶


















In [47]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_colu_2014', 'BF-09_Houet_colu_2012']"
contig="3L"



     



















In [48]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [49]:



     


window=1000



     



















In [50]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [51]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_coluzzii_Chr3L.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_coluzzii_Chr3L.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_coluzzii_Chr3L.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_coluzzii_Chr3L.tsv




























COLUZZII¶
3R¶


















In [52]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_colu_2014', 'BF-09_Houet_colu_2012']"
contig="3R"



     



















In [53]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [54]:



     


window=1000



     



















In [55]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [56]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_coluzzii_Chr3R.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_coluzzii_Chr3R.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_coluzzii_Chr3R.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_coluzzii_Chr3R.tsv




























COLUZZII¶
X¶


















In [57]:



     


sample_query="cohort_admin2_year == ['BF-09_Houet_colu_2014', 'BF-09_Houet_colu_2012']"
contig="X"



     



















In [58]:



     


ag3.plot_h12_calibration(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    


















    
    







  











    
    
































In [59]:



     


window=5000



     



















In [60]:



     


ag3.plot_h12_gwss(
    contig=contig,
    window_size=window,
    sample_sets=sample_sets,
    sample_query=sample_query,
)



     






















    
    















    
    


















    
    







  











    
    
































In [61]:



     


# Parameters
output_file_path = os.path.join(results_dir, 'H12counts_coluzzii_ChrX.tsv')
outliers_file_path = os.path.join(results_dir, 'H12counts_95qoutliers_coluzzii_ChrX.tsv')


# Perform H12 calculation
h12counts = ag3.h12_gwss(
    contig=contig,
    sample_sets=sample_sets,
    sample_query=sample_query,
    window_size=window,
)

# Initialize an empty list to store the results
results = []

# Process h12counts based on its actual type and structure
if isinstance(h12counts, tuple):
    window_positions, h12_values = h12counts
    # Convert window positions to integers
    window_positions = window_positions.astype(int)
else:
    print("Unexpected h12counts format")

# Calculate start and end positions based on the window size
start_positions = window_positions - (window // 2)
end_positions = window_positions + (window // 2)

# Save the results to the output file in tabular format with headers and contig name
with open(output_file_path, 'w') as file:
    # Write the headers
    file.write("chr\tstart\tend\tmid\tH12\n")
    # Write data
    for start, end, mid, h12 in zip(start_positions, end_positions, window_positions, h12_values):
        file.write(f"{contig}\t{start}\t{end}\t{mid}\t{h12}\n")

print(f"Results saved to {output_file_path}")

# Load the results into a DataFrame
df = pd.read_csv(output_file_path, delimiter='\t')

# Compute the 95th percentile for H12 values
threshold = np.percentile(df['H12'], 95)

# Filter the DataFrame to get the top 5% H12 values
outliers_df = df[df['H12'] >= threshold]

# Save the outliers to a new file
outliers_df.to_csv(outliers_file_path, sep='\t', index=False)

print(f"Outliers saved to {outliers_file_path}")



     






















    
    






Results saved to drive/MyDrive/malariagen_data_cache/H12counts_coluzzii_ChrX.tsv
Outliers saved to drive/MyDrive/malariagen_data_cache/H12counts_95qoutliers_coluzzii_ChrX.tsv
























In [69]:



     


%cd /content/drive/MyDrive/Colab\ Notebooks/



     






















    
    






/content/drive/MyDrive/Colab Notebooks
























In [68]:



     


# Step 2: Convert the notebook to HTML
import os

# Define the filename
notebook_filename = "drive/MyDrive/Colab\ Notebooks/H12stats_Burkina.ipynb"
html_filename = notebook_filename.replace(".ipynb", ".html")

# Save the notebook
!jupyter nbconvert --to html "{notebook_filename}"

# Step 3: Download the HTML file
from google.colab import files

# Download the HTML file
files.download(html_filename)



     






















    
    






[NbConvertApp] WARNING | pattern 'drive/MyDrive/Colab\\ Notebooks/H12stats_Burkina.ipynb' matched no files
This application is used to convert notebook files (*.ipynb)
        to various other formats.

        WARNING: THE COMMANDLINE INTERFACE MAY CHANGE IN FUTURE RELEASES.

Options
=======
The options below are convenience aliases to configurable class-options,
as listed in the "Equivalent to" description-line of the aliases.
To see all configurable class-options for some <cmd>, use:
    <cmd> --help-all

--debug
    set log level to logging.DEBUG (maximize logging output)
    Equivalent to: [--Application.log_level=10]
--show-config
    Show the application's configuration (human-readable format)
    Equivalent to: [--Application.show_config=True]
--show-config-json
    Show the application's configuration (json format)
    Equivalent to: [--Application.show_config_json=True]
--generate-config
    generate default config file
    Equivalent to: [--JupyterApp.generate_config=True]
-y
    Answer yes to any questions instead of prompting.
    Equivalent to: [--JupyterApp.answer_yes=True]
--execute
    Execute the notebook prior to export.
    Equivalent to: [--ExecutePreprocessor.enabled=True]
--allow-errors
    Continue notebook execution even if one of the cells throws an error and include the error message in the cell output (the default behaviour is to abort conversion). This flag is only relevant if '--execute' was specified, too.
    Equivalent to: [--ExecutePreprocessor.allow_errors=True]
--stdin
    read a single notebook file from stdin. Write the resulting notebook with default basename 'notebook.*'
    Equivalent to: [--NbConvertApp.from_stdin=True]
--stdout
    Write notebook output to stdout instead of files.
    Equivalent to: [--NbConvertApp.writer_class=StdoutWriter]
--inplace
    Run nbconvert in place, overwriting the existing notebook (only
            relevant when converting to notebook format)
    Equivalent to: [--NbConvertApp.use_output_suffix=False --NbConvertApp.export_format=notebook --FilesWriter.build_directory=]
--clear-output
    Clear output of current file and save in place,
            overwriting the existing notebook.
    Equivalent to: [--NbConvertApp.use_output_suffix=False --NbConvertApp.export_format=notebook --FilesWriter.build_directory= --ClearOutputPreprocessor.enabled=True]
--no-prompt
    Exclude input and output prompts from converted document.
    Equivalent to: [--TemplateExporter.exclude_input_prompt=True --TemplateExporter.exclude_output_prompt=True]
--no-input
    Exclude input cells and output prompts from converted document.
            This mode is ideal for generating code-free reports.
    Equivalent to: [--TemplateExporter.exclude_output_prompt=True --TemplateExporter.exclude_input=True --TemplateExporter.exclude_input_prompt=True]
--allow-chromium-download
    Whether to allow downloading chromium if no suitable version is found on the system.
    Equivalent to: [--WebPDFExporter.allow_chromium_download=True]
--disable-chromium-sandbox
    Disable chromium security sandbox when converting to PDF..
    Equivalent to: [--WebPDFExporter.disable_sandbox=True]
--show-input
    Shows code input. This flag is only useful for dejavu users.
    Equivalent to: [--TemplateExporter.exclude_input=False]
--embed-images
    Embed the images as base64 dataurls in the output. This flag is only useful for the HTML/WebPDF/Slides exports.
    Equivalent to: [--HTMLExporter.embed_images=True]
--sanitize-html
    Whether the HTML in Markdown cells and cell outputs should be sanitized..
    Equivalent to: [--HTMLExporter.sanitize_html=True]
--log-level=<Enum>
    Set the log level by value or name.
    Choices: any of [0, 10, 20, 30, 40, 50, 'DEBUG', 'INFO', 'WARN', 'ERROR', 'CRITICAL']
    Default: 30
    Equivalent to: [--Application.log_level]
--config=<Unicode>
    Full path of a config file.
    Default: ''
    Equivalent to: [--JupyterApp.config_file]
--to=<Unicode>
    The export format to be used, either one of the built-in formats
            ['asciidoc', 'custom', 'html', 'latex', 'markdown', 'notebook', 'pdf', 'python', 'rst', 'script', 'slides', 'webpdf']
            or a dotted object name that represents the import path for an
            ``Exporter`` class
    Default: ''
    Equivalent to: [--NbConvertApp.export_format]
--template=<Unicode>
    Name of the template to use
    Default: ''
    Equivalent to: [--TemplateExporter.template_name]
--template-file=<Unicode>
    Name of the template file to use
    Default: None
    Equivalent to: [--TemplateExporter.template_file]
--theme=<Unicode>
    Template specific theme(e.g. the name of a JupyterLab CSS theme distributed
    as prebuilt extension for the lab template)
    Default: 'light'
    Equivalent to: [--HTMLExporter.theme]
--sanitize_html=<Bool>
    Whether the HTML in Markdown cells and cell outputs should be sanitized.This
    should be set to True by nbviewer or similar tools.
    Default: False
    Equivalent to: [--HTMLExporter.sanitize_html]
--writer=<DottedObjectName>
    Writer class used to write the
                                        results of the conversion
    Default: 'FilesWriter'
    Equivalent to: [--NbConvertApp.writer_class]
--post=<DottedOrNone>
    PostProcessor class used to write the
                                        results of the conversion
    Default: ''
    Equivalent to: [--NbConvertApp.postprocessor_class]
--output=<Unicode>
    overwrite base name use for output files.
                can only be used when converting one notebook at a time.
    Default: ''
    Equivalent to: [--NbConvertApp.output_base]
--output-dir=<Unicode>
    Directory to write output(s) to. Defaults
                                  to output to the directory of each notebook. To recover
                                  previous default behaviour (outputting to the current
                                  working directory) use . as the flag value.
    Default: ''
    Equivalent to: [--FilesWriter.build_directory]
--reveal-prefix=<Unicode>
    The URL prefix for reveal.js (version 3.x).
            This defaults to the reveal CDN, but can be any url pointing to a copy
            of reveal.js.
            For speaker notes to work, this must be a relative path to a local
            copy of reveal.js: e.g., "reveal.js".
            If a relative path is given, it must be a subdirectory of the
            current directory (from which the server is run).
            See the usage documentation
            (https://nbconvert.readthedocs.io/en/latest/usage.html#reveal-js-html-slideshow)
            for more details.
    Default: ''
    Equivalent to: [--SlidesExporter.reveal_url_prefix]
--nbformat=<Enum>
    The nbformat version to write.
            Use this to downgrade notebooks.
    Choices: any of [1, 2, 3, 4]
    Default: 4
    Equivalent to: [--NotebookExporter.nbformat_version]

Examples
--------

    The simplest way to use nbconvert is

            > jupyter nbconvert mynotebook.ipynb --to html

            Options include ['asciidoc', 'custom', 'html', 'latex', 'markdown', 'notebook', 'pdf', 'python', 'rst', 'script', 'slides', 'webpdf'].

            > jupyter nbconvert --to latex mynotebook.ipynb

            Both HTML and LaTeX support multiple output templates. LaTeX includes
            'base', 'article' and 'report'.  HTML includes 'basic', 'lab' and
            'classic'. You can specify the flavor of the format used.

            > jupyter nbconvert --to html --template lab mynotebook.ipynb

            You can also pipe the output to stdout, rather than a file

            > jupyter nbconvert mynotebook.ipynb --stdout

            PDF is generated via latex

            > jupyter nbconvert mynotebook.ipynb --to pdf

            You can get (and serve) a Reveal.js-powered slideshow

            > jupyter nbconvert myslides.ipynb --to slides --post serve

            Multiple notebooks can be given at the command line in a couple of
            different ways:

            > jupyter nbconvert notebook*.ipynb
            > jupyter nbconvert notebook1.ipynb notebook2.ipynb

            or you can specify the notebooks list in a config file, containing::

                c.NbConvertApp.notebooks = ["my_notebook.ipynb"]

            > jupyter nbconvert --config mycfg.py

To see all available configurables, use `--help-all`.











    
    






---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-68-bbca913fce8c> in <cell line: 15>()
     13 
     14 # Download the HTML file
---> 15 files.download(html_filename)

/usr/local/lib/python3.10/dist-packages/google/colab/files.py in download(filename)
    223   if not _os.path.exists(filename):
    224     msg = 'Cannot find file: {}'.format(filename)
--> 225     raise FileNotFoundError(msg)  # pylint: disable=undefined-variable
    226 
    227   comm_manager = _IPython.get_ipython().kernel.comm_manager

FileNotFoundError: Cannot find file: drive/MyDrive/Colab\ Notebooks/H12stats_Burkina.html